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Pure Monte Carlo Method: a Third Way for Plasma Simulation 
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We bring a totally new concept for plasma simulation, other than the conventional two ways: 
Fluid/Kinetic Continuum (FKC) method and Particle-in-Cell (PIC) method. This method is based 
on Pure Monte Carlo (PMC), but far beyond traditional treatments. PMC solves all the equations 
(kinetic, fluid, fleld) and treats all the procedures (collisions, others) in the system via MC method. 
As shown in two paradigms, many advantages have found. It has shown the capability to be the 
third importance approach for plasma simulation or even completely substitute the other two in the 
future. It's also suitable for many unsolved problems, then bring plasma simulation to a new era. 

PACS numbers: 52.65.-y, 52.65.Pp 



Introduction-The original of Monte Carlo (MC) 
method for computations can back to hundreds of years 
ago, e.g., the Buffon's needle experiment. With the emer- 
gence of the modern computer in 1940s and pioneer works 
by John von Neumann, Stanislaw Ulam and Nicholas 
Metropolis, this method brings amounts of practical ap- 
plications. The famous and classical [Mctropolisl953] pa- 
per is said Q mainly contributed by the hero of plasma 
physics, M. N. Rosenbluth. However, up to now, MC 
for physics is mainly on statistical physics, and only in 
merely special cases for plasma physics. 

Some people had known that, in fact, the PIC method 
can be seen a type of MC method. [l| provides a uni- 
fied MC interpretation of particle simulations, which then 
can use the MC theory to estimate the errors and develop 
many types of sampling approaches, e.g., simple MC, im- 
portance sampling, the S-i method and more. 

While, at most, the conventionally PIC only treats 
the kinetic equation using MC viewpoint and the field 
equations are still solved by Deterministic Discretiza- 
tions (DD, finite differences and so on). FKC is solved 
totally by DD. Generally, each approach has its advan- 
tages and disadvantages. The main disadvantages of DD 
are the numerical unstability and dissipation and not 
easy to treat complex geometries and complex procedures 
which cannot be described easily by differential equations 
(ODEs, PDEs). To avoid old disadvantages, a way out 
is to develop new approaches. 

We find, indeed, we can have a third way for plasma 
simulation, i.e., we use totally stochastic statistical 
(Monte Carlo) approach. We need only a change of con- 
cept. In Copenhagen interpretation of quantum mechan- 
ics, we have known that the intrinsic of the world is uncer- 
tainty and all things are probabilities, only the statistic 
is determined. If we can accept this, we can easily ac- 
cept the philosophy of this PMC letter. The remaining 
thing is how we do that, or exactly, how we solve differ- 
ential equations (in fluid and kinetic) via MC. A branch 
of mathematic called Stochastic Differential Equations 
(SDE) has provided us some basic tools. 

At below, we will firstly describe the PMC method and 



then treat a fluid edge pedestal transport model and the 
kinetic ID electrostatic problem as paradigms. 

PMC Method-ln plasma simulation, MC has been 
used to many complex procedures, such as collisions and 
impurities transport, which has been widely known. So, 
here, we need not talk too much about that, and all the 
old MC methods in plasma simulation can be inherited 
by PMC method naturally. The new steps of PMC are 
to tell how we solve the fluid and field equations. 

Using MC to solve PDEs can back to two pivotal pa- 
pers, [J| and IJ], the former uses probabilistic interpreta- 
tions for linear elliptic and parabolic equations, the latter 
gives MC methods for solving integro-difFcrcntial equa- 
tions occur in various branches of the natural sciences. 

Feynman-Kac formula-Many types of plasma equa- 
tions (e.g., transport equations, Vlasov equation) have 
the form: 
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Feynman-Kac formula [6[[l0| tells us, Eq.((T]) is equivalent 
to SDE when V^p — Q {V^p^Q case can also find). 



dX = ^iiX, t)dt + a{X, t)dWt, 



(2) 



where Wt is Wiener process (Brownian motion), which 
is described by normal (Gaussian) distribution. Eqs.([T]) 
and ^ are easily general to high dimensions or more 
variables, e.g., (x, v,t). The above descriptions are also 
shown by Ito's lemma Q, a classical work in SDE. 

Now, we solve the convection-diffusion equation as ex- 
ample to show how to use this formula, which is stan- 
dard steps and well known by SDE people but may not 
familiar by plasma community. At this case, a = /i and 
D = are constants. For /(x, 0) = /o(a;) , wc start 

with a set samples ■Ci i ' ■ ■ from foix). A sample 
for every next step is 



t+dt 



^* -f adt + V2Ddir]„ i = l,--- ,N. (3) 



where rji is normal distribution with zero mean and unit 
variance. For a ~ 0, we get pure diffusion, and for D = 
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0; pure convection. In fact, integral form explicit exact 
solution is 

f{x,t) = J Gt{x - y)fo{y - at)dy = Gt * fo{x - at), 

(4) 



MC (|)(x,y)=x^+y^ 
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If we use MC to integral (jU directly, and use N' samples 
for particles and M points in space, the total cost will 
be 0{N' X M) . While using dS]), the total cost is 0{N). 
We write Q here for benchmark. 

Using ([3]), a result is show in FIG[TJ We can see that 
the MC result reproduces the exact solution very well, 
when we using larger N and M, the result can be even 
better. Here, M is not for calculation but just for recon- 
structing the grids for /. 
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FIG. 1: (Color online) Solve convection-diffusion equation us- 
ing (|3|, as example to show PMC method. 

The above is a summary of SDE for our usage, which 
is a key fundament of our PMC method, the applications 
to plasma examples will be shown at below. 

Poisson equation-Wc may meet Poisson equation fre- 
quently in plasma physics, for example, in the elec- 
trostatic kinetic case or the electromagnetic gyrokinetic 
case. The equation is as the form 



(x) ==-p(x), xe D. 



(5) 



If p = 0, gives Laplace equation, which can be found 
solved by MC in many places, e.g., 1^. A good in- 



troduction of MC methods for Poisson equation can 
find in [11 . The method can also take from an ex- 
tended Feynman-Kac formula. For Dirichlet boundary, 
(/)(x) = /(x), x G dD. Standard solution for position xo 
using SDE notation is 
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piWt)dt 



+ E[f{Wr,^)], (6) 



where, tqd = inf{t : Wt £ dD} is the first-passage time 
and Wrgo is the first-passage location on the boundary, 
E is average. Using floating random walk, the results 
for 2D are show in FIGH] The MC method here has 
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FIG. 2: (Color online) MC for 2D Poisson equation, example. 



another advantage that we can obtain the solution at a 
few points directly instead of calculating the whole range, 
and even if there are steep gradients or irregular bound- 
ary. Some updates of MC for Poisson equation may also 
find in recent literatures, such as for Neumann boundary, 
improving efficiencies, or reducing errors. 

Other equations-There are many other types of PDEs 
in plasma physics. Scalar conservation law equation 
du/dt + dF{u)/dx = Oje.g., Burgers' equation, can find 
be solved generally in [Oj. Navier-Stokes (vector) equa- 
tion for turbulence can find in For ODE, we often 
need just the standard MC integral. Some more MC 
methods and details for electromagnetic problems are in- 
troduced in [l^ . One can also refer [llj for some general 
descriptions. 

The types of equations mentioned at above have con- 
tented many of the main plasma physics equations. 

Application Paradigms~We use the above methods 
to solve two practical plasma problems. 

Simple fluid edge pedestal transport model— The sim- 
plest model for calculating pedestal width can find in 
[lit or [2^ , the inner range of SOL (Scrape-Off Layer, for 
divertor or limiter) is mainly slowly perpendicular trans- 
port, and the outer range is mainly fast parallel (to the 
magnetic field line, mainly on toroidal f direction) trans- 
port and the particles travel at most L {2'kR or 2TiqR) 
then hit the target with the thermal velocity, where R is 
the major radius, q is edge safety factor, and the ther- 
mal velocity is equal to acoustic velocity at edge from 
sheath calculations. 

Wc can write the inner range transport equation as 
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, _ (7) 
For simplification, we only use the density diffusion equa- 
tion and treat the coefficient D as constant, and also ig- 
nore the source term S. For outer range, we find quickly 
that we cannot easily write a ID equation combine to 
the inner range equation. Even though we write down 
the outer range equation, wc will quickly meet the singu- 
larity problem at the last closed (magnetic) fiux surface 
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(LCFS, separatrix). A special treatment for transport 
barriers can find in [l9| . which rewrites the transport 
equations to new form. However, using PMC, we can 
avoid the above problems automatically. 

The inner equation ([7]) can be solved using ([2]). The 
outer range needs an extra MC treatment: if a particle 
comes out to the outer range, we use the same equation 
as the inner range for perpendicular transport but in the 
same time let it vanish random in time [0, T) if it is still 
in the outer range, with T = L/cg. 

Quickly, we can find when using ([3]) to treat (O, the 
drift velocity a = —D/r — >■ oo for the MC particles at 
r — )• 0. This problem is also met similarly in PIC Q when 
using polar coordinate, also in large scale PIC simulation 
as in GTC The conventional method to treat this 

problem is making a hollow and ignoring a small r ^ 
region. While, one can find in [l^ other ways of MC im- 
plement transport equations and treat the r — > prob- 
lem (e.g., L'Hospital's rule limp^^df /{pdp) = d"^ j jdf?), 
which means that this won't be a severe problem in PMC. 

However, in fact, we will find MC can be a natural 
method to solve this geometry and complex boundary 
problems. The only reason why people using different 
coordinates is that we can use the symmetry of the sys- 
tem to make us treat (e.g., analytical calculations) or 
understand problem easier. People care but the great na- 
ture never cares the coordinates. The physics law won't 
change in any coordinates. So, wc can rewrite ([7]) from 
(r, 9) coordinate to Cartesian coordinate (x, y) and use 
\J -V ip' = r = rt) as boundary. MC is suitable for any 
complex boundaries, and the circle boundary is just one 
of them. Then, in use of PMC, people can use Carte- 
sian coordinate always (in fact, as the circle symmetry 
case in FIGHJ we have already given an example), which 
will simplify many old complex coordinates (notable: the 
magnetic surface coordinate) or boundaries problems a 
lot. A 2D implement of PMC for this edge pedestal trans- 
port model is shown in FIG13] In tokamak experiment. 



the system is described by kinetic equation 
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FIG. 3: (Color online) PMC for (Tokamak) edge pedestal 
transport. 

one can see lots of figures as FIGlU i-C-, the above simple 
model can model the experiment at least qualitatively. 

Kinetic ID electrostatic (ESID) problem-FoT ESID, 



dt ' ^ dx ^ m dv 



(8) 



where, q is the charge, m is mass, E ~ — V0 is the elec- 
trical field, C(/) is collision term. The field equation is 
dl]), with p = J {fi — fe)dv the charge density. All the 
variables has normalized in units of e(j>/Te, w~e^, Xdc- For 
ID, ([5]) can rewrite to a more simple ODE form 

dE/dx = p. (9) 

Firstly, we can easily find ^ has the form of ([1]), gives 

dxi = Vidt, dvi = {qE/m)dt, (10) 



which is totally the same as in PIC. Eg. ([10)) tells us that 
more dimensions won't bring significant new time cost 
(note: higher dimensions we may need more particles for 
accuracy, which will bring some new time cost), while 
will if we use DD. If the collision term has the form 
C(/) = af + I3df /dv+'-fd^ f /dv^ , the kinetic equation is 
still nothing else but ([T]). We see here, for PIC simulation 
of collision, we need an extra MC step, while in PMC 
method, we can write this step to equation of motion 
(ITUl) directly. 

Here, we assume ions immobile for simplify and treat 
the coUisionless (Vlasov) case, i.e., C(/) ~ , because 
that which is easy for benchmark. 

For f[x,v,t) , X and v can be anywhere. For E(yx) or 
(/)(a;) , the x can also be anywhere. But, ([8|) and (|9]) are 
coupled. We should find out a way to connect the posi- 
tion X in / and E. A naturally thought is calculating the 
E{x) directly using MC for every /-particle one by one, 
but which is too much time consumed, especially when 
/-particles [Nj) and 0-particles [N^) are very large. An 
adjust way is using grids and then interpolating, which 
can reduce the time cost by an order of 0{N^). We find 
this is just what PIC does, i.e., we re-deduce all the key 
steps of PIC method by PMC! The interpolating method 
is easy and has checked by PIC, however, one may also 
find other ways to implement this PMC step if can make 
sure they work. 

Using ODE ([9]) as field equation, for MC, which is just 
an integral from one side to another side, and we use ran- 
dom particles to do this, the errors of each step will cancel 
and average is zero. While, for DD (e.g., Euler, Runge- 
Kutta or else), the integral is summation and errors will 
accumulate. Using Poisson equation, the MC method is 
provided at before, an extra step is needed to calculate E 
from (j) by calculating gradient from the MC ^-particles. 
Periodic boundary for the MC particles of / is easy to 
implement, which is the same as PIC. We comment here, 
for (j), the periodic boundary means (/)(0) = 0(i), but 
when we see E, which is not E{0) = E{L) but average 
< E{x) >= (Q has known this). In fact, this implies 
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TABLE I: Comparison of FKC, PIC and PMC for kinetic plasma simulation. 





FKC 


PIC 


PMC 




f ,riTii~i ni 1 1 1 TTi 


Pfirticlp -1- pnntiniiiiTTi 


Onlv nflrtiplp 


Equations 


Only PDEs 


New method for eqn. of motion, 
Field eqn. brown from FKC 


Eqn. of motion brown from PIC, 
New method for field equation 


Good 


Accurate 


Efficient 


Easy for complex problems. Effi- 
cient for parallel. Easy for coding 


Bad 


Numerical unstability and dissipation 


Noise 


Crude, with error oc l/\/N 


Philosophy 


More mathematical 


Half math, half first principle 


More first principle 



that, we have added two charged slabs (not real 'parti- 
cles', because this is ID, or called 'markers' in PIC) at 
the boundary to cancel the extra E field! Few ESID PIC 
researchers have seen this. Thinking in PMC way will 
bring us many new insights or re-thoughts. 

For this ESID problem, we can find, the conventional 
PIC treatment and new PMC treatment can be very sim- 
ilar. The implements can be close, but the concepts are 
totally different. PIC is fixed and nonadjustable. PMC 
provides more possibilities and can bring us new ways. 

Using PMC to solving ^ and ([9]), a result is show in 
FIG. m Note that, due to noises, the phase space plot- 
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FIG. 4: (Color online) Compare with PMC and PIC for 2- 
stream instability, phase space plotting. 

ting will change for both PIC and PMC, even though 
under the same parameters (e.g., keep source code un- 
changed), especially in nonlinear stages. The above fig- 
ure is a select from typical runs. The main feature of 
phase space holes is clearly in both PMC and PIC re- 
sults. PMC is slightly slower than PIC in this case. 

Summary and Comments-The above two 
paradigms have shown many new advantages. In 
fact, PMC is more powerful when problems are com- 
plicated. For example, many large scaled codes are 
possible be completely rewritten via this new method, 
an example is BOUT ([Ilj, gyrofiuid). For large 
scaled gyro-kinetic, this method is also possible. Some 
complicated paradigms are beyond this short letter. 

A comparison of FKC, PIC and PMC is hsted in 
TABEL. U which can also shows that why we can call 
PMC be a third method for plasma simulation. 

At present, many tools of PMC are brown from other 



fields. But the concept change is the key important thing. 
The unified viewpoint by combining everything in one 
provides us a totally new picture. In this new angle, we 
can see difficult old problems naturally and problems be- 
come simple. We can not only do new problems but also 
provide new insights or more efficient methods for old 
problems. For kinetic problems, PMC can be seen as ad- 
vanced PIC; for fluid problems and complex procedures, 
PMC is new. Whatever, people may just concern that 
if an approach is useful and works well. PMC matches 
this. PMC is not a single fixed approach as PIC, but a 
combine of varies of old and new MC methods, flexible. 

As far as the author (HSX) knows, Max-Planck- 
Institut fiir Plasmaphysik is developing MC methods for 
Stellar ator edge transport physics, while it seems [l^l 
has already or almost implemented a fully PMC ES2D 
code though they didn't know this concept when they 
did that. New applications for resistive MHD tearing 
mode, Hasegawa-Mima equation, Grad-Shafranov equa- 
tion and many other problems are on the way by the 
author or others. This letter has given/built the frame- 
work of PMC, which is enough for many practical appli- 
cations. But, there may are still many new intelligent 
works need do for some special types of cases. Some un- 
expected problems may also meet when using to more 
complicated examples. The author does not know that 
yet. If also implemented well for the above untested sev- 
eral examples, this new approach is possible completely 
substitute the conventional two methods, FKC and PIC, 
and also suitable for solving many old unsolved or hard 
to solve problems, then bring plasma simulation to a new 
era. 

Acknowledgement-This work is also inspired by Zhi- 
chen FENG (IFTS-ZJU), who once (2010) proposed a 
new method for PIC, though which performances not well 
in practical application yet. 
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